/*
    Copyright 2007-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Implementation and explicit instantiations of the intensity_stage
/// class. \ingroup mcrx

// $Id$

#include "config.h"
#include "mcrx-stage.h"
#include "mcrx-stage-impl.h"
#include "create_grid.h"
#include "arepo_grid.h"
#include "arepo_emission.h"

using namespace std;

template <template <typename> class grid_type> 
bool mcrx::intensity_stage<grid_type>::check_stage_state () const
{
  bool state = false;
  try {
    this->output_file_->pHDU().readKey("MCRX" + stage_ID(), state);
  }
  catch (...) {}
  return state;
}

template <template <typename> class grid_type> 
void mcrx::intensity_stage<grid_type>::set_stage_state (bool state)
{
  this->output_file_->pHDU().addKey("MCRX" + stage_ID(), state, "");
}

template <template <typename> class grid_type> 
void mcrx::intensity_stage<grid_type>::setup_objects ()
{
  // read original and intensity wavelength vectors
  CCfits::ExtHDU& lambda_hdu = open_HDU (*this->output_file_, "LAMBDA");
  array_1 lambda_orig;
  read(lambda_hdu.column("lambda"), lambda_orig);
  read(lambda_hdu.column("lambda_intensity"), lambda_);
  int n_lambda_int;
  lambda_hdu.readKey("NLAMBDA_INTENSITY", n_lambda_int);
  lambda_.resizeAndPreserve(n_lambda_int);

  // Build emission, asking it to subsample to new grid.
  cout << "Setting up emission" << endl;
  CCfits::ExtHDU& pdata_hdu = open_HDU(*this->input_file_, "PARTICLEDATA");
  this->emi_.reset(new T_emission(pdata_hdu, 
				  this->p_.defined("use_reference_for_emission") &&
				  this->p_.getValue("use_reference_for_emission", bool()),
				  lambda_orig, lambda_));

  // Set up dust model
  cout << "Setting up dust model" << endl;

  // Read dust grains
  typename T_dust_model::T_scatterer_vector 
    sv(read_dust_grains<typename T_dust_model::T_scatterer_vector>(this->p_, this->m_.units));
  this->model_.reset(new T_dust_model (sv));
  this->model_->set_wavelength(lambda_);

  // Now build grid
  cout << "Building grid" << endl;
  
  // The density generator is used to generate the densities of the
  // dust grains from the gas and metal densities in the simulation.  
  uniform_density_generator dg(this->p_);

  // by definition, we track intensities in this stage
  full_sed_grid<grid_type>* dummy=0;
  this->g_ = create_grid(*this->input_file_, dg, lambda_.size(), dummy);

  // The intensity stage doesn't use any cameras, so we just generate
  // an empty emergence object
  this->eme_.reset(new T_emergence (vector<std::pair<T_float, T_float> >(),
			      1,1,1));

  // find reference wavelength
  vector<T_float> templambda(lambda_.begin(), lambda_.end());
  reflambda_= 
    lower_bound (templambda.begin(), templambda.end(),
		 this->p_.getValue("reference_wavelength", T_float ())) -
    templambda.begin();
  assert(reflambda_>=0);
  assert(reflambda_<lambda_.size());

  cout << "Reference wavelength: " << lambda_ (reflambda_) << " ("
       << reflambda_<< ")\n";
}


template <template <typename> class grid_type> 
void mcrx::intensity_stage<grid_type>::load_file ()
{
  // first set up general stuff
  setup_objects();

  // Only thing to load is grid intensity
  this->g_->load_intensity(open_HDU(*this->output_file_, "INTENSITY"));
}


template <template <typename> class grid_type> 
void mcrx::intensity_stage<grid_type>::load_dump (binifstream& dump_file)
{
  // first set up general stuff
  setup_objects();

  // Only thing to load is grid intensity
  this->g_->load_dump(dump_file);
}


template <template <typename> class grid_type> 
void mcrx::intensity_stage<grid_type>::save_file ()
{
  // Get normalization and save cell intensities
  const T_float normalization = 1./this->n_rays_;
  this->g_->write_intensity(*this->output_file_, 
			    "INTENSITY","DEPOSITION",
			    normalization, *this->model_,
			    this->p_.defined("compress_images")&&
			    this->p_.getValue("compress_images", bool ()) );

  if(this->model_->absorbed().size()>0)
    cout << "Luminosity removed from ray tracing (absorbed): " 
	 << integrate_quantity(this->model_->absorbed(),this->lambda_)*normalization;
}


template <template <typename> class grid_type> 
void mcrx::intensity_stage<grid_type>::save_dump (binofstream& dump_file) const
{
  this->g_->write_dump(dump_file);
}


template <template <typename> class grid_type> 
bool mcrx::intensity_stage<grid_type>::shoot ()
{
  return this->shoot_scatter (true);
}

template <template <typename> class grid_type> 
void mcrx::intensity_stage<grid_type>::operator() ()
{
  // this is here so the function is instantiated
  this->run_stage();
}


// explicit instantiations
template class mcrx::intensity_stage<mcrx::adaptive_grid>;

#ifdef WITH_AREPO
template class mcrx::intensity_stage<mcrx::arepo_grid>;
#endif
